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1.  Introduction 


One  of  the  fundamental  problems  facing  discrimination  and  location  of  earthquakes  and 
explosions  is  creating  accurate  structure  models  for  the  crust  and  mantle.  The  routine  practice  of 
locating  seismic  events  based  on  a  one-dimensional  velocity  model  inevitably  introduces  bias 
into  locations.  Using  a  high-resolution  three-dimensional  (3D)  seismic  velocity  model 
significantly  improves  seismic  event  location  accuracy  (Smith  and  Ekstrom,  1996;  Firbas,  2000; 
Steck  et  al.,  2000)  and  thus  helps  to  satisfy  the  goal  of  nuclear  explosion  monitoring.  In  addition 
to  seismic  event  locations,  seismic  wave  amplitudes  are  also  important  in  discriminating  between 
earthquakes  and  explosions.  For  example,  the  Fg/Pg  ratio  has  been  shown  to  be  an  effective 
discriminant  (Jiao  et  al.,  2003).  As  seismic  waves  travel  through  an  anelastic  and  heterogeneous 
medium,  their  amplitudes  will  be  attenuated  and  knowledge  of  the  attenuation  structure  is  vital  to 
correct  the  distortion.  Here  we  developed  high-resolution  models  of  the  seismic  velocity  and 
attenuation  structure  of  the  Sichuan- Yunnan  region  (97°-108°E  and  21°-35°N)  that  is  located  in 
southwest  China  using  seismic  catalog  and  waveform  data. 

The  Sichuan- Yunnan  region  lies  in  an  active  transition  zone  between  the  Yangtze  platform  to  the 
east  and  the  Tibetan  plateau  to  the  west  and  is  generally  believed  to  have  had  its  origin  in  the 
collision  between  the  Indian  plate  and  the  Eurasian  plate  about  45  Ma  ago.  This  continental 
collision  led  to  active  tectonic  deformation  on  a  larger-scale  and  created  a  high  level  of 
seismicity  (Figure  1).  There  are  seven  major  active  seismic  zones  in  this  region  including  the 
Fongmen  Shan  belt,  the  Xianshuihe  seismic  belt,  the  Anninghe  seismic  belt,  the  Xiaojiang 
seismic  belt,  the  Red  River  seismic  belt,  the  Fancang-Gengma  seismic  belt,  and  the  Tengchong- 
Fongling  seismic  belt  (Chan  et  al.,  2001).  Most,  but  not  all,  of  the  earthquakes  in  these  belts  are 
associated  with  fault  zones  that  are  identified  at  the  surface  (Figure  1).  Obtaining  a  high- 
resolution  seismic  velocity  model  for  this  region  and  precise  earthquake  locations  will  help 
delineate  fault  zones  at  depth  more  clearly  and  facilitate  their  association  with  surface  fault 
traces.  Several  magnitude  7  earthquakes  occurred  in  this  region  in  the  last  twenty  years.  Recent 
examples  are  the  1995  Mengnian  earthquake  (M=7.4)  and  the  1996  Fijiang  earthquake  (M=7.0). 
The  auxiliary  IMS  station  KMI  is  located  in  this  region  and  China’s  nuclear  test  site  Fop  Nor  is 
located  about  1000  km  to  the  northwest. 
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Figure  1.  Relocated  earthquakes  in  the  Yunnan  (red  dots)  and  Sichuan  (black  dots)  region  from 
this  study.  White  and  green  triangles  are  Yunnan  and  Sichuan  provincial  seismic  stations. 

Black  lines  are  mapped  fault  traces  on  the  surface.  FI:  Xiaojiang  Fault;  F2:  Red  River  Fault; 
F3:  Wuliangshan  Fault;  F4:  Lancangjiang  Fault;  F5:  Nandinghe  Fault;  F6:  Nujiang  Fault;  F7: 
Qiaohou-Weishan  Fault;  F8:  Tonghai-Chuxiong  Fault;  F9:  Anninghe  Fault;  F10:  Zemuhe 
Fault;  FI  1 :  Longmenshan  Fault;  F12:  Xianshuihe  Fault;  F13:  Jinshajiang  Fault. 

Most  of  the  previous  seismic  tomography  studies  that  cover  the  Sichuan- Yunnan  region  are  on  a 
global  scale.  Examples  include  the  CUB  1.0  (Shapiro  and  Ritzwoller,  2002)  and  the  Science 
Applications  International  Corporation  (SAIC)  global  l°xl°  models  that  were  constructed  from 
group  and  phase  velocity  dispersion  measurements  of  surface  waves,  and  the  CRUST  5.1 
(Mooney,  1998)  and  CRUST  2.0  (Laske  et  ah,  2001)  models  from  seismic  refraction  data.  There 
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are  also  some  regional  models  in  China  such  as  the  Pn  model  of  Hearn  and  Ni  (2001),  the  Pn  and 
Sn  velocity  model  of  Ritzwoller  et  al.  (2002),  P-wave  models  of  Liu  et  al.  (1990),  Xu  et  al. 

(2002),  and  Sun  et  al.  (2004),  and  surface-wave  models  of  Wu  et  al.  (1997)  and  Lebedev  and 
Nolet  (2003).  The  grid  spacing  ranges  from  l°xl°  to  5°x5°  and  are  too  coarse  to  characterize  the 
detailed  structure  of  the  Sichuan- Yunnan  region.  The  most  recent  work  by  Huang  et  al.  (2003) 
and  Liang  et  al.  (2004)  used  a  smaller  grid  spacing  of  0.5°  to  calculate  a  Pn  model  for  this  region 
that  shows  the  Pn  velocity  of  uppennost  mantle  varies  from  7.7~8.3  km/s.  Regional  3D  velocity 
models  are  also  available  for  the  Yunnan-Sichuan  region  from  several  studies  (He  et  al.,  2005; 
Huang  et  al,  2002;  Liu  et  al,  2005;  Wang  et  al.,  2003).  Huang  et  al.  (2002)  determined  a  3-D  Vp 
model  for  the  Sichuan-Yunnan  region  (22°-32°N,  98°-105°E)  with  a  horizontal  grid  spacing  of 
50  to  80  km  and  a  vertical  grid  spacing  of  10-25  km  using  1315  local  and  regional  earthquakes  of 
magnitude  >  2.5  and  172  stations.  S  times  were  only  used  to  locate  earthquakes.  Wang  et  al. 
(2003)  used  P  and  S  arrival  times  from  4625  local  and  regional  earthquakes  recorded  at  174 
stations  to  determine  Vp  and  Vs  models  of  southwestern  China  (21°  34°N,  97°-105°E).  The  grid 
spacing  in  horizontal  directions  is  0.5°  and  the  vertical  spacing  ranges  from  10  km  to  20  km. 
Different  from  Huang  et  al.  (2002)  and  Wang  et  al.  (2003),  Liu  et  al.  (2005)  used  arrival  times 
from  602  local  and  regional  earthquakes  and  985  teleseismic  events  to  detennine  a  3D  Vp  model 
beneath  the  north-south  tectonic  belt  between  Tibet  and  Eastern  China.  The  grid  spacing  in  the 
horizontal  direction  is  0.5°  and  the  vertical  spacing  is  20-35  km  except  for  a  finer  grid  interval  of 
3  km  near  the  surface.  In  addition  to  3D  seismic  tomography  studies  using  body  waves,  He  et  al. 
(2005)  used  short-period  (1-18  s)  surface  wave  data  recorded  by  23  stations  in  the  Yunnan  region 
to  determine  a  3D  S-wave  velocity  structure  of  the  middle  and  upper  crust. 

Our  recent  development  and  applications  of  double-difference  (DD)  tomography  showed  that  the 
new  method  has  the  ability  to  characterize  the  detailed  velocity  structure  and  improve  seismic 
event  locations  in  a  superior  manner  compared  to  conventional  seismic  tomography  (Zhang  and 
Thurber,  2003,  2005,  2006;  Zhang  et  al.,  2004;  Thurber  et  al.,  2004).  DD  tomography  is  a 
generalization  of  DD  location  (Waldhauser  and  Ellsworth,  2000);  it  simultaneously  solves  for  the 
3D  velocity  structure  and  seismic  event  locations.  DD  tomography  uses  a  combination  of 
absolute  and  more  accurate  differential  arrival  times  and  hierarchically  detennines  the  velocity 
structure  from  larger  scale  to  smaller  scale.  Recent  application  of  DD  location  in  this  region 
(Yang  et  al.,  2003;  2005)  showed  that  the  relocated  events  are  more  concentrated  near  the  fault 
zones  using  just  the  catalog  differential  data.  From  their  location  results,  we  are  confident  that 
even  using  the  catalog  data  alone  will  significantly  improve  the  velocity  model,  though  the  best 
results  will  be  achieved  with  wavefonn  alignment  methods. 

To  our  knowledge,  there  has  been  no  attenuation  tomography  study  published  for  this  region. 
There  are  some  surface  wave  (Lg)  and  Pn  wave  attenuation  studies  in  the  Tibetan  plateau  that  is 
located  to  the  west  of  our  targeted  region  (Xie  et  al.,  2002;  Xie,  2003).  They  found  low  Lg  Q 
values  (at  1  Hz)  of  typically  lower  than  200  in  the  entire  plateau  that  are  consistent  with  high 
temperature  and/or  fluids  in  the  Tibetan  crust  (Xie  et  al.,  2002).  The  frequency  independent  P 
wave  Q  values  imply  that  a  melt-bearing  uppermost  mantle  might  lie  underneath  Tibet  (Xie, 
2003). 

Because  of  the  high  seismicity  (more  than  1000  earthquakes  with  magnitudes  greater  than  2.5 
per  year)  and  severe  earthquake  hazards,  the  Sichuan  and  Yunnan  Seismological  Bureaus  have 
been  operating  dense  local  networks  of  stations  since  the  1970s  (Figure  1).  Most  of  the  stations 
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were  upgraded  to  short-period  and  broadband  (up  to  30  s)  digital  stations  after  1999.  We 
acquired  P  and  S  arrival  times  recorded  by  local  network  stations  for  the  entire  Sichuan- Yunnan 
region  (Figure  1).  We  also  acquired  seismic  waveforms  for  the  Yunnan  region  in  the  period  of 
1999  to  2004.  Using  these  catalog  picks  and  waveforms,  we  determined  3D  Vp  and  Vs  models 
for  the  Yunnan- Sichuan  region  and  3D  Qp  and  Qs  models  for  the  Yunnan  region.  Our  research 
project  consists  of  four  components: 

(1)  using  waveform  alignment  methods  (waveform  cross-correlation  (WCC)  and  bispectrum  (BS) 
analysis)  to  obtain  more  accurate  differential  arrival  times  for  the  Yunnan  region; 

(2)  using  regional  scale  adaptive-grid  DD  tomography  to  obtain  detailed  P-  and  S-wave  velocity 
models  of  the  Sichuan- Yunnan  region; 

(3)  developing  an  adaptive-grid  “triple-difference”  seismic  attenuation  method  and  applying  it  to 
determine  the  detailed  seismic  attenuation  structure  for  both  Qp  and  Qs  for  the  Yunnan 
region; 

(4)  assembling  a  ground  truth  database. 

In  the  following  sections,  we  will  give  a  detailed  description  of  each  component  including 
methods  and  results. 


2,  Data  collection  for  the  Yunnan-Sichuan  region 

Yunnan  region.  We  collected  3-component  waveforms  recorded  by  26  Yunnan  provincial 
stations  for  -5000  events  in  the  period  of  1999  to  2004.  It  was  time  consuming  and  laborious  to 
process  these  waveforms  to  be  used  for  WCC  analysis.  First  we  converted  all  the  waveforms 
stored  in  the  original  China  data  fonnat  into  standard  SAC  format  required  by  the  BCSEIS 
package  (Du  et  al.,  2004).  Then  we  matched  these  waveforms  with  the  catalog  picks  we  already 
collected  and  filled  the  corresponding  SAC  headers.  For  the  waveforms  that  did  not  have 
matched  catalog  picks,  we  first  used  the  wavelet- AIC  picker  (Zhang  et  al.,  2003)  to 
automatically  pick  first  P  and  S  arrivals.  The  picker  was  originally  designed  to  only  pick  first  P 
arrivals.  We  adapted  it  to  pick  first  S  arrivals  on  horizontal  components.  We  first  calculated  the 
first  S  arrivals  based  on  the  catalog  locations  and  an  approximate  3D  Vs  model.  Then  we  started 
the  time  window  for  the  wavelet- AIC  picker  5  seconds  before  the  estimated  S  arrivals  and 
compared  the  picks  selected  for  the  two  horizontal  components.  If  the  two  selected  picks  are 
reasonably  close  to  each  other  (e.g.  within  -1  second),  then  we  chose  the  mean  of  the  two  picks 
as  the  first  S  arrivals.  For  quality  control,  we  visually  checked  all  the  waveforms  and  corrected 
and  removed  those  outliers.  In  this  way,  we  can  guarantee  the  quality  of  P  and  S  picks.  After  this 
process,  we  obtained  -35000  P  and  -30000  S  picks  for  -5600  events  recorded  by  26  stations 
(Figure  1). 

Sichuan  region.  We  assembled  and  analyzed  catalog  picks  for  various  local  networks  in  the 
Sichuan  region.  In  total,  there  are  -  41400  P  and  -38700  S  picks  corresponding  to  -5000 
earthquakes  and  -50  stations  (Figure  1).  Compared  to  the  Yunnan  region,  there  were  not  many 
waveforms  available.  We  first  analyzed  a  total  of  2651  high-quality  waveforms  for  267  events  at 
16  stations.  We  further  analyzed  an  additional  -80  events  that  do  not  have  catalog  locations  and 
origin  times.  We  estimated  initial  locations  and  origin  times  for  these  events  using  a  grid  search 
method  NonLinLoc  developed  by  Anthony  Lomax.  The  absolute  P  and  S  arrivals  extracted  from 
these  waveforms  along  with  other  catalog  picks  are  used  for  seismic  tomography. 
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3.  Waveform  cross-correlation  using  the  bispectrum  (BS)  method,  BCSEIS 

Researchers  using  cross-correlation  (CC)  to  align  waveforms  often  choose  those  time  delay 
estimates  with  CC  coefficients  above  a  specified  threshold.  For  example,  Schaff  et  ah,  (2002) 
only  select  those  time  delays  with  CC  values  larger  than  0.70  and  mean  coherences  above  0.70. 
The  selection  of  an  optimum  threshold  value  for  waveform  cross-correlation  is  important  but 
difficult.  If  it  is  set  too  high,  then  only  a  limited  number  of  very  accurate  differential  time  data 
are  available  for  further  analysis.  If  the  threshold  value  is  set  too  low,  then  many  unreliable 
differential  time  estimates  are  used  which  will  negatively  affect  the  relocation  and  tomography 
results. 

The  BCSEIS  algorithm  that  we  adopted  for  measuring  wavefonn  CC  times  works  in  the  third- 
order  spectral  domain  and  can  suppress  correlated  Gaussian  or  low-skewness  noise  sources 
(Nikias  and  Raghuveer,  1987;  Nikias  and  Pan,  1988).  BCSEIS  employs  a  third-order  spectral 
method  to  calculate  two  additional  time  delay  estimates  with  both  the  raw  (unfiltered)  and  band¬ 
pass  filtered  waveforms,  and  uses  them  to  verify  (select  or  reject)  the  one  computed  with  the  CC 
technique  using  the  filtered  waveforms  (Du  et  ah,  2004).  Thus  this  BS  verification  process  can 
reject  unreliable  CC  time  delay  estimates  and  also  can  accept  additional  CC  time  delays  even  if 
their  associated  CC  coefficients  are  smaller  than  a  nominal  threshold  value  if  they  pass  the  BS 
verification  procedure  (Du  et  ah,  2004). 

We  have  calculated  differential  times  for  -38000  waveforms  for  -5600  events  recorded  by  26 
stations  in  Yunnan  Province  using  BCSEIS  and  obtained  -130000  P-  and  -33000  S-wave 
differential  travel  times  given  estimated  origin  times.  Figure  2  shows  the  distribution  of  the 
calculated  P-  and  S-wave  differential  times.  From  the  figure  we  can  see  that  the  majority  of 
differential  times  are  less  than  1  second.  Such  measurements  are  associated  with  the  event  pairs 
with  a  separation  of  a  few  kilometers.  The  other  measurements  with  larger  differential  times 
correspond  to  those  event  pairs  that  could  be  separated  by  more  than  50  km  and  can  also  have 
similar  P  and  S  waveforms.  We  found  that  the  calculated  differential  times  with  different 
filtering  are  slightly  different  for  some  waveform  pairs.  To  assure  that  the  results  are  not  affected 
by  filtering,  we  tested  several  different  filters.  We  found  that  the  default  Butterworth  filter  with  3 
poles  and  2  passes  in  BCSEIS  changes  the  waveform  shape  around  the  first  arrivals  significantly. 
After  testing  different  filtering,  we  decided  to  choose  a  Butterworth  filter  with  2  poles  and  1  pass 
in  our  correlation  computation.  We  also  tested  different  window  sizes  to  find  a  suitable  one  for 
our  analysis. 

Figures  3  and  4  show  the  improvement  of  CC  differential  times  over  catalog  picks  for  those 
similar  events.  In  Figure  3,  a  set  of  similar  waveforms  observed  at  the  station  BST  are  aligned 
with  the  catalog  P  arrivals  and  CC  differential  times,  separately.  The  stacked  waveforms  are 
shown  in  the  bottom  panels.  Note  that  the  waveforms  are  aligned  relatively  well  using  catalog 
picks,  verifying  the  quality  of  the  catalog  picks  (Figure  3).  For  comparison,  we  can  see  that  the 
quality  of  the  CC  alignment  is  much  better  than  that  of  the  catalog  pick  aligmnent.  As  shown  in 
the  figures  by  vertical  lines,  the  time  window  we  choose  here  is  30  samples  before  and  97 
samples  after  the  preliminary  P  arrival  picks.  Similarly,  Figure  4  shows  wavefonn  alignments 
with  catalog  S  arrivals  and  CC  differential  times  for  another  set  of  similar  wavefonns  recorded  at 
the  station  HQT.  Due  to  the  greater  scattering  (errors)  in  S  arrival  picking,  the  improvement  of 
the  CC  aligmnent  for  S  phases  is  even  more  significant.  Similarly,  we  also  obtained  a  total  of 
-4100  P  and  S  differential  times  from  323  events  recorded  by  17  stations  in  the  Sichuan  region. 
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number  or  measurement 


Figure  2.  Distribution  of  the  calculated  P-  and  S-wave  differential  times  using  the  waveform 
cross-correlation  package  BCSEIS  for  the  Yunnan  region. 


Figure  3.  Waveform  alignment  according  to  (left)  catalog  P  picks  and  (right)  waveform  cross¬ 
correlation  P  times. 
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Figure  4.  Waveform  alignment  according  to  (left)  catalog  S  picks  and  (right)  waveform  cross¬ 
correlation  S  times. 
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4.  Regional-scale  adaptive-grid  double-difference  tomography  method 

Adaptive  tomography  method  based  on  tetrahedral  diagrams 

The  two  most  commonly  used  schemes  to  represent  the  Earth  in  seismic  tomography  are  the 
regular  constant-velocity  cell  approach  (Aki  and  Lee,  1976)  and  the  regular  3D  grid  approach 
(Thurber,  1983).  The  advantage  of  the  regular  constant-velocity  cell  approach  is  its  simplicity, 
but  because  the  velocity  is  constant  inside  a  cell  and  discontinuous  from  one  cell  to  another,  cells 
are  not  ideal  for  approximating  a  smoothly  varying  structure,  and  finding  accurate  ray  paths 
connecting  events  and  stations  can  be  difficult.  For  the  regular  3D  grid  approach,  the  velocity 
varies  continuously  among  nodes  in  all  directions  using  linear  B-spline  interpolation  (Thurber, 
1983)  or  cubic  B-spline  interpolation  (Michelini  and  McEvilly,  1991).  No  matter  which  approach 
is  taken,  however,  the  ray  distribution  typically  is  highly  uneven  due  to  inadequate  acquisition 
geometry,  uneven  distribution  of  seismic  sources,  missing  data,  and  ray  bending.  Some  grid 
points  or  cells  may  have  few  or  even  no  rays  sampling  them,  while  other  cells  may  have  very 
dense  rays  sampling  them.  The  regular  grid  or  cell  spacing  restriction  makes  it  difficult  to  adapt 
the  model  to  the  uneven  distribution  of  ray  paths.  The  mismatch  between  the  ray  distribution  and 
the  grid  or  cells  used  results  in  instability  of  the  inversion. 

The  goal  of  matching  the  ray  distribution  and  the  inversion  grid  or  cells  in  seismic  tomography 
naturally  leads  to  the  irregular  grid  or  cells  distribution.  Sambridge  and  Gudmundsson  (1998) 
and  Bohm  et  al.  (2000)  proposed  irregular  cell  approaches  on  the  basis  of  tetrahedral  and 
Voronoi  diagrams,  respectively.  In  the  former  study,  they  did  not  explicitly  explain  how  to 
optimize  the  grid.  In  the  latter  study,  the  cells  are  adjusted  based  on  the  model  null  space  energy 
(Vesnaver,  1996)  and/or  velocity  gradient.  For  whole  Earth  tomography,  Sambridge  and  Faletic 
(2003)  recently  proposed  a  data-driven  tetrahedral  cell  adaptive  scheme  based  on  the  maximum 
spatial  gradients  in  seismic  velocity  perturbation  across  each  tetrahedron  face.  The 
tetrahedral/Voronoi  diagram  methods  allow  a  more  flexible  representation  of  the  model.  For 
example,  cells  of  widely  varying  sizes  with  complex  distributions  are  easily  implemented,  and  it 
is  more  convenient  to  build  parameterizations  containing  particular  interfaces,  on  which  the 
nodes  can  be  distributed. 


Figure  5.  Tetrahedral  diagram  for  a  set 
of  8  nodes,  defining  a  total  of  6 
tetrahedra. 


Figure  6.  Set  of  tetrahedral  before  (left)  and 
after  (right)  adding  a  node.  Only  the 
tetrahedron  within  which  the  new  node  is 
inserted  is  affected  by  the  process. 
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We  adopt  the  tetrahedral  parameterization  approach  and  associated  interpolation  schemes  for  our 
adaptive-grid  DD  tomography  algorithm.  Ideally,  the  distribution  of  the  inversion  grid  should 
match  with  the  data  distribution:  in  regions  where  there  are  more  rays,  the  inversion  grid  nodes 
should  be  more  finely  distributed.  The  inversion  problem  will  then  be  better  conditioned.  Here 
we  use  the  derivative  weight  sum  (DWS)  values  at  the  grid  nodes  to  indicate  the  ray  density 
(Thurber,  1983).  We  start  the  inversion  from  a  regular  inversion  grid,  equivalent  to  that  of  the 
simul2000  algorithm  (Thurber  and  Eberhart-Phillips,  1999).  We  randomly  perturb  the  starting 
regular  inversion  grid  by  a  very  small  amount  (so  that  the  nodes  are  not  located  on  the  same 
plane)  in  order  to  construct  the  tetrahedral  or  Voronoi  diagram  around  these  nodes  using  the 
Qhull  algorithm  (Barber  et  al.,  1996).  Figure  5  shows  the  tetrahedral  diagram  for  8  starting 
regular  nodes.  We  also  construct  a  regular  computational  grid  that  remains  fixed  during  the 
inversion,  which  can  be  the  same  as  the  starting  regular  inversion  grid  or  finer,  following  the 
strategy  of  Kissling  et  al.  (2001).  On  a  local  scale,  we  trace  rays  between  events  and  stations 
based  on  the  current  regular  velocity  grid  using  a  Cartesian  pseudo-bending  method  (Um  and 
Thurber,  1987).  The  rays  between  all  the  event  and  station  pairs  are  saved  for  later  use  in 
defining  the  adaptive  grid.  Using  these  rays,  we  find  the  partial  derivatives  of  the  travel  times 
with  respect  to  the  model  slowness  parameters  on  the  initial  (slightly  irregular)  inversion  grid.  In 
the  process,  we  calculate  the  DWS  values  on  the  inversion  grid  nodes.  We  then  automatically  go 
through  each  tetrahedron  to  check  the  sum  of  the  DWS  values  at  its  4  nodes.  If  the  sum  is  larger 
than  a  predefined  threshold  threshold  1,  an  additional  node  will  be  inserted  at  the  center  of  this 
tetrahedron.  The  velocity  value  at  this  node  is  interpolated  by  the  linear  interpolation  method 
from  the  4  inversion  grid  nodes  or  by  the  natural  neighbor  interpolation  method  from  its  N 
natural  neighbors  (Sambridge  et  al.,  1995),  depending  on  which  interpolation  method  has  been 
adopted.  Figure  6  shows  three  tetrahedra  constructed  from  5  nodes.  One  additional  node  is 
inserted  into  the  middle  of  the  bottom  tetrahedron  creating  4  new  tetrahedra  (Figure  6).  The 
advantage  of  refining  one  tetrahedron  by  adding  a  node  inside  it  is  to  avoid  affecting  other 
neighboring  tetrahedra.  At  the  same  time,  we  also  remove  all  the  inversion  grid  nodes  having  a 
DWS  value  smaller  than  the  threshold  threshold2.  Generally  thresholdl  is  30  times  larger  than 
threshold2.  Note  that  we  keep  all  the  irregular  grid  nodes  that  were  originally  the  boundary  grid 
nodes  of  the  starting  regular  inversion  grid  to  guarantee  all  the  regular  computational  grid  nodes 
are  located  inside  the  tetrahedral  or  Voronoi  diagram. 

After  this  refining  process,  we  obtain  a  new  irregular  inversion  grid.  This  irregular  inversion  grid 
cannot  be  guaranteed  to  be  data-adaptive  without  further  quality  control.  We  construct  a  revised 
tetrahedral  or  Voronoi  diagram  around  the  new  irregular  inversion  grid  and  recalculate  the  DWS 
values  on  all  the  nodes  using  the  previously  saved  ray  information.  After  this,  we  examine  all  the 
irregular  nodes  to  check  if  the  DWS  value  on  each  one  is  greater  than  the  predefined  value 
threshold2.  The  irregular  grid  nodes  not  satisfying  the  requirement  will  be  removed.  This  new 
irregular  grid  is  the  inversion  grid  to  be  used  for  the  current  iteration  of  simultaneous  inversion. 
Finally  a  new  tetrahedral  diagram  is  constructed  and  the  partial  derivatives  of  the  travel  times 
with  respect  to  the  new  set  of  inversion  grid  nodes  are  determined  for  the  construction  of  the 
seismic  tomography  equations.  After  each  simultaneous  inversion,  the  velocity  values  on  the 
irregular  inversion  grid  nodes  and  the  regular  computational  grid  nodes  are  updated.  For 
subsequent  simultaneous  inversions,  the  inversion  grid  is  again  adaptively  updated  to  better 
match  with  the  data  distribution,  following  the  same  procedure. 
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“Sphere-in-a-box  ’’finite-difference  method 

In  the  “sphere-in-a-box”  version  of  DD  tomography,  the  curvature  of  the  Earth  is  explicitly  taken 
into  account.  We  note  that  the  earth  flattening  transformation  approach  is  only  valid  for  ID 
velocity  models  and  not  for  3D  velocity  models  so  that  approach  is  not  a  viable  one.  Following 
Flanagan  et  al.  (2007),  we  solve  this  problem  by  parameterizing  a  spherical  surface  inside  a 
Cartesian  volume  of  grid  nodes.  Figure  7a  shows  the  case  of  putting  the  entire  Earth  into  a  cube 
so  the  curvature  of  the  Earth  can  be  taken  into  account.  This  scheme  can  be  used  for  global 
seismic  tomography.  For  regional  seismic  tomography,  we  are  only  interested  in  a  part  of  the 
Earth  and  we  can  construct  a  rectangular  box  covering  the  portion  of  interest  (Figure  7b  and  c). 
The  coordinate  center  is  placed  at  the  surface  of  the  Earth,  positive  X  and  Y-axes  point  to  the 
direction  of  North  and  West,  and  positive  Z-axis  points  downward.  The  grid  nodes  above  the 
Earth’s  surface  (air  nodes)  are  given  the  velocity  for  P  waves  traveling  in  air  (Figure  7c).  As  a 
result,  all  the  rays  travel  inside  the  Earth.  Here  we  should  emphasize  that  this  version  is  still 
based  on  a  Cartesian  coordinate  system  and  not  on  a  true  spherical  system. 


Figure  7.  Grid  setup  for  regional 
scale  double-difference 
tomography,  (a)  The  whole 
earth  is  inserted  into  a  cubic  box. 
(b)  A  rectangular  box  covers  the 
region  of  interest,  (c)  The 
representation  of  (b)  in  the  X-Z 
plane.  The  surface  of  the  Earth 
is  shown  as  a  thick  line.  The 
regular  grid  is  used  for  finite- 
difference  ray  tracing  algorithm. 


Finite-difference  ray  tracing  algorithms  require  a  uniformly  spaced  velocity  (or  slowness)  grid 
(Podvin  and  Lecomte,  1991).  We  interpolate  the  velocity  values  on  the  regular  uniform  grid 
nodes  from  non-uniform  inversion  grid  nodes  through  trilinear  interpolation.  First  we  treat  each 
station  as  a  source  and  calculate  travel  times  to  all  velocity  nodes  in  the  volume.  The  travel  time 
from  a  station  to  each  earthquake  is  interpolated  from  its  8  neighboring  nodes  through  trilinear 
interpolation.  The  ray  path  from  the  earthquake  to  the  station  is  found  iteratively  with  increments 
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opposite  to  the  travel  time  gradient.  The  corresponding  partial  derivatives  of  travel  times  with 
respect  to  the  slowness  models  are  calculated  by  dividing  the  ray  paths  into  small  segments 
(Thurber,  1983). 

Key  modifications  of  the  tomography  code 

During  the  project,  we  made  several  key  modifications  to  make  the  code  applicable  to  large  data 
set  on  regional  scales.  First  we  modified  the  adaptive-mesh  tomography  code  from  local  scale  to 
regional  scale  by  including  the  “sphere-in-a-box”  finite-difference  travel  time  calculation  method. 
Previously  the  code  had  difficulty  saving  ray  path  information  to  adapt  the  inversion  mesh  due  to 
the  limitation  of  memory  space  (Zhang  and  Thurber,  2005).  We  overcame  this  challenge  by 
modifying  the  code  to  save  ray  density  information  on  an  intennediate  grid  and  then  project  it 
back  to  the  adaptive  inversion  mesh. 

We  also  modified  the  code  by  finding  a  way  to  apply  the  smoothing  constraint  to  the  inversion 
mesh.  Previously,  we  set  up  an  intermediate  regular  grid  and  then  projected  the  smoothing 
weights  on  it  to  the  inversion  mesh  (Zhang  and  Thurber,  2005).  This  could  deteriorate  the 
advantage  of  using  an  adaptive  mesh  to  resolve  the  fine  scale  structure  because  the  spatial  scale 
of  the  structure  is  limited  by  the  size  of  the  intennediate  grid.  We  now  have  found  a  way  to 
efficiently  search  the  natural  neighbors  for  each  inversion  mesh  node  and  are  able  to  directly 
apply  the  first-order  smoothing  constraint. 

The  adaptive-mesh  DD  tomography  code  is  further  modified  by  also  adding  new  nodes  at  the 
midpoints  of  tetrahedral  edges  when  a  tetrahedron  is  sampled  by  many  rays,  as  determined  by  a 
threshold  level.  The  previous  strategy  was  to  only  add  new  nodes  to  the  center  of  one  tetrahedron. 
The  new  strategy  will  be  helpful  for  more  uniformly  distributing  the  inversion  mesh  nodes. 


5.  Seismic  attenuation  tomography 

Seismic  attenuation  can  provide  relatively  independent  and  complementary  information  on 
subsurface  structure  (Sanders,  1993;  Eberhart-Phillips  et  ah,  1995;  Romero  et  ah,  1997). 
Examples  of  attenuation  tomography  applied  to  fault  zone  settings  can  be  found  in  Lees  and 
Lindley  (1994),  for  Loma  Prieta,  and  Rietbrock  (2001),  for  Kobe. 

Our  approach  will  follow  that  of  Rietbrock  (2001).  Briefly,  the  set  of  spectra  for  a  given  event  at 
each  of  the  observing  stations  j  was  fit  to  an  Q-square  source  model  with  a  single  comer 
frequency,  using 

ln(Aj(f))  =  ln(Qoj)  -  ln(  1  +  (f/fc)2)  -Tift,*  (1) 

where  Aj  is  the  spectral  amplitude  at  frequency  f  for  the  event  observed  at  station  j,  Q0j  is  the 
long-period  spectral  plateau  level,  fc  is  the  corner  frequency  for  the  event,  and  tj*  is  the  whole 
path  attenuation  operator  for  the  event  path  to  station  j  (Rietbrock,  2001): 

t*.  =  f  dr  (2) 

U  JpalV1J(r)Qij(r) 
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Equation  (1)  is  fit  for  the  range  of  frequencies  with  adequate  signal-to-noise  ratio.  The  set  of  all 
t*  values  is  then  inverted  for  the  3D  Q  structure,  using  the  3D  seismic  velocity  model  and 
associated  event  locations  to  determine  the  ray  paths  for  the  t*  integral.  In  fact,  once  the  event 
locations  and  the  3D  velocity  model  are  determined,  the  relationship  between  the  t*  values  and 
the  reciprocal  of  the  Q  structure  (1/Q)  is  linear. 

Spectral-ratios  between  P  and  S  arrivals  at  one  station  or  between  P  arrivals  at  two  stations  for  a 
single  event  have  been  used  to  derive  differential  t*  values  that  are  then  inverted  for  relative  Q 
variations  (Roth  and  Wiens,  1999).  This  is  essentially  a  “double-difference”  attenuation 
tomography  method  (difference  between  observed  and  calculated  differences  in  t*),  which 
removes  some  of  the  systematic  components  from  the  equation.  This  can  be  taken  one  step 
further  to  produce  a  “triple-difference”  attenuation  tomography  method,  similar  to  the  approach 
of  Blakeslee  et  al.  (1989).  Expressing  the  amplitude  spectrum  from  source  i  to  station  m  as 

Ahn(f)=  Ai(f>0im^im)  W*>  Glm«  ^(^0  (3) 


where  A;  is  the  radiated  source  spectrum,  Im  is  the  instrument  response,  and  G,m  is  the  geometric 
spreading,  then  for  a  pair  of  nearby  events  i  and  j  observed  at  a  pair  of  nearby  stations  1  and  m, 
we  form  the  ratio  of  ratios 


Anffl/A,mffl 


=  C(r) 


exp(-7tf(t*  +  t*m )) 
exp(-7tf(t*m  + 1  *j )) 


(4) 


where  C(r)  is  a  ratio  of  ratios  of  the  geometrical  spreading  terms  (assumed  to  be  independent  of 
frequency),  with  the  instrument  responses  and  source  spectral  amplitudes  canceling  out 
(assuming  that  the  radiated  source  spectra  from  a  given  event  at  two  nearby  stations  are 
approximately  equal).  Taking  the  natural  log,  we  arrive  at 


Vf)/Ajm(f)J 


=  In 


(5) 


which  is  a  linear  function  of  frequency  and  can  be  fit  to  determine  the  observed  t*  difference 
term.  Then  the  difference  between  the  observed  and  calculated  t*  difference  terms  can  be  related 
to  perturbations  to  the  Q  model  via 


[fij’l  -til)-(t*n  -Cl]*  "[(tj,  -t*,)-(t*m  -Of'  =  J 


dr 


dr 


V„AQ,  JV,AQ, 


J - — - J 

J  V  AO  j 


dr 


V  AO  J  V  AO 

jm  Xjm  im  ^mit 


(6) 


The  left-hand  side  can  be  thought  of  as  a  "triple  difference."  In  principle,  the  difference  data  will 
be  more  accurate,  and  we  can  combine  t*  residuals,  differential  t*  residuals  (i.e.  site-response 
corrected),  and  "triple  difference"  data  to  derive  a  more  accurate  Q  model  in  a  manner  analogous 
to  DD  seismic  velocity  tomography  (discussed  above).  We  developed  this  refined  attenuation 
modeling  approach  as  part  of  this  project.  The  adaptive-grid  approach  proposed  in  the  above 
section  is  also  adapted  in  the  attenuation  tomography  in  the  same  way  as  we  do  for  the  seismic 
tomography. 
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Measurement  of  the  absolute  t*  Values 

We  fit  spectral  amplitude  data  to  equation  (1)  within  the  range  of  frequencies  with  adequate 
signal-to-noise  ratio  to  calculate  t*  values.  We  use  the  multitaper  method  of  spectral  analysis 
(Park  et  ah,  1987)  to  calculate  amplitude  spectra  from  the  windowed  vertical  waveforms  around 
the  P  arrivals.  We  also  estimate  the  noise  spectra  from  the  seismograms  in  a  window  right  before 
the  signal  window. 

In  the  traditional  single-taper  analysis,  portions  of  the  waveforms  of  interest  are  excluded  from 
analysis  as  a  trade-off  for  reducing  the  spectral  leakage.  In  the  multitaper  approach,  a  family  of 
statistically  independent  spectral  estimates  is  computed  from  a  signal  using  a  set  of  orthogonal 
tapers  that  are  referred  to  as  discrete  prolate  spheriodal  sequences  (DPSS).  Averaging  over  this 
ensemble  of  spectra  yields  an  estimate  with  much  lower  variance  than  that  from  single-taper 
methods.  Stable  spectral  estimation  is  of  importance  for  a  robust  measurement  of  t*  values. 

Figure  8  shows  samples  of  seismograms  recorded  at  two  stations  HQT  and  EYA  for  two  events 
with  magnitude  3.3  and  4.0,  respectively.  The  signal  windows  used  for  amplitude  spectrum 
calculations  are  indicated  by  vertical  lines  in  the  figure.  The  corresponding  multitaper  spectra 
estimated  from  these  seismograms  are  shown  in  Figure  9.  In  each  panel  of  the  figure,  the 
calculated  amplitude  spectra  are  shown  as  solid  lines,  while  the  noise  spectra  are  represented  by 
dotted  lines.  The  fits  of  spectral  amplitude  data  to  equation  (1)  are  demonstrated  with  dashed 
lines.  We  obtain  a  set  of  CL,fc  and  t*  values  from  the  fits  as  shown  in  each  panel  of  the  figure 
(Figure  9).  Using  this  method,  we  calculated  -10200  P-wave  t*  values  from  -3000  events  and 
-8860  S-wave  t*  values  from  -2300  events  in  the  Yunnan  region. 
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Figure  8.  Waveforms  at  two  nearby  stations  EYA  and  HQT  from  two  nearby  events. 
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Figure  9.  Spectral  fitting  based  on  Equation  1  to  estimate  t*  values.  Blue  solid  line:  original 
spectra.  Red  dashed  line:  fitting  spectra.  Blue  dots:  noise  spectra. 


Measurement  of  differential  t*  values  using  two  events  and  two  stations 

The  above  approach  for  measuring  absolute  t*  values  may  potentially  suffer  limitations  when 
there  is  significant  site  response  at  some  stations.  We  have  proposed  to  use  two  events  and  two 
stations  to  remove  source  signature  and  site  response  from  the  spectral  amplitude  using  equations 
(4)  and  (5). 

Figure  10  shows  the  calculated  spectral  ratios  from  the  above  two  stations  and  two  events.  We  fit 
spectral  ratios  to  the  linear  equation  (5)  and  get  the  differential  t*  value  of  0.0066.  The 
differential  t*  value  calculated  directly  from  the  measured  absolute  t*  values  is  0.0070.  These 
two  values  are  close;  however,  the  measured  differential  t*  values  from  spectral  ratios  are  free 
from  station  effects  and  are  not  affected  by  the  source  model  assumption.  These  differential  t* 
values  can  be  used  to  solve  for  the  Q  model  in  a  way  similar  to  DD  tomography.  We  note  that 
the  multitaper  method  also  improves  the  robustness  of  spectral  ratios  significantly. 


13 


As  can  be  seen  in  Figure  10,  the  estimated  differential  t*  values  may  be  affected  significantly  by 
the  fitting  frequency  range.  We  developed  an  automatic  process  to  calculate  the  estimated 
differential  t*  values  in  which  the  fitting  frequency  range  is  chosen  based  on  signal-to-noise  ratio 
(SNR)  and  fitting  error.  The  optimal  frequency  range  is  selected  when  the  SNR  is  larger  than  3 
and  the  standard  deviation  of  the  fitting  line  to  the  spectra  is  smallest.  We  measured  -120,000 
differential  t*  values  for  P  waves  using  this  automatic  process  for  -1550  events  on  23  stations  of 
the  Yunnan  seismic  network.  Figure  1 1  shows  the  comparison  of  differential  P-wave  t* 
measurements  with  respect  to  those  calculated  directly  using  the  Q-square  source  model.  It 
indicates  that  when  the  frequency  range  is  larger,  the  discrepancy  is  smaller.  For  example,  for  the 
fitting  frequency  range  greater  than  37  Hz,  there  is  a  clear  linear  trend  between  the  two 
measurements,  indicating  both  absolute  and  differential  t*  measurements  are  reliable.  Using  the 
same  procedure,  we  also  calculated  -122,500  differential  S-wave  t*  values  for  1455  events  on  24 
stations.  We  compared  these  differential  S-wave  t*  values  to  those  calculated  directly  from  the 
absolute  t*  values  using  the  Q-square  source  model.  A  similar  conclusion  to  that  for  P  waves  can 
also  be  drawn. 


Figure  10.  Differential  t*  measurements  using  spectral  ratios  for  two  events  and  two  stations 
shown  in  Figure  9. 
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Figure  11.  Comparison  of  measured  differential  t*  values  and  those  calculated  directly  from 
absolute  t*  values  for  P  waves. 


6.  Earthquake  locations  and  3D  velocity  models  in  the  Yunnan  region 

For  the  Yunnan  region,  we  obtained  -96000  first  P  arrivals  and  -32000  first  S  arrivals  for  -5600 
events  and  26  stations  (Figure  12).  From  them,  we  constructed  -289000  P  and  -254000  S 
differential  times  for  event  pairs  within  40  km.  A  maximum  of  20  neighbors  is  allowed  when 
forming  links  among  event  pairs.  The  average  number  of  links  (or  common  observations)  per 
event  pair  is  8  and  the  average  inter-event  distance  is  1 1  km.  Along  with  -96000  P  and  -32000  S 
cross-correlation  times,  we  used  the  double-difference  tomography  method  to  simultaneously 
determine  earthquake  locations  and  the  Vp  and  Vs  models. 

We  took  a  multi-step  strategy  to  update  event  locations  and  the  velocity  models.  First  we  started 
from  a  regular  inversion  grid  in  latitude,  longitude  and  depth  and  an  estimate  ID  model.  The  grid 
spacing  in  latitude  and  longitude  is  0.5°  and  it  is  5  km  in  depth  from  0  km  to  80  km.  The  initial 
RMS  residuals  for  catalog  times  and  CC  times  are  3.059  s  and  5.471  s,  respectively.  We  adopted 
a  hierarchical  weighting  system  for  different  data  types  (Zhang  and  Thurber,  2003,  2006).  We 
first  applied  10  times  higher  weighting  to  the  absolute  times  to  image  the  larger  scale  structure 
and  move  earthquakes  closer  to  their  locations.  Then  we  applied  higher  weighting  to  the 
differential  times  (including  catalog  differential  and  CC  times)  to  further  refine  the  structure  near 
the  source  region  and  improve  earthquake  locations.  The  final  absolute  and  weighted  RMS 
residuals  for  catalog  times  are  569  ms  and  356  ms.  For  comparison,  the  final  absolute  and 
weighted  RMS  residuals  for  CC  times  are  69  ms  and  39  ms.  Both  smoothing  (weight=30)  and 
damping  are  applied  to  stabilize  the  inversion.  A  grid  interval  of  5  km  is  used  for  finite- 
difference  travel  time  calculation. 
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Figure  12.  Event  and  station  distribution  for  the  Yunnan  region.  FI:  Xiaojiang  Fault;  F2:  Red 
River  Fault;  F3:  Wuliangshan  Fault;  F4:  Lancangjiang  Fault;  F5:  Nandinghe  Fault;  F6: 
Nujiang  Fault;  F7:  Qiaohou-Weishan  Fault;  F8:  Tonghai-Chuxiong  Fault;  F9:  Anninghe  Fault; 
F10:  Zemuhe  Fault;  Rl:  Tengchong-Longling  seismic  belt;  R2:  Lancang-Gengma  seismic  belt; 
R3:  Lijiang  seismic  belt;  R4:  Yongshen-Ninglang-Muli-Juilong  seismic  belt;  and  R5:  Tonghai- 
Shiping  seismic  belt. 

We  then  started  from  the  earthquake  locations  and  velocity  models  inverted  from  the  previous 
step  and  used  only  the  more  accurate  CC  times  to  further  update  the  velocity  models  using  the 
adaptive  grid  approach.  Figure  13  shows  the  grid  nodes  at  depths  of  7.5  to  12.5  km  after  the 
inversion.  It  is  clear  to  see  that  the  density  of  grid  nodes  follows  the  earthquake  distribution 
pattern  shown  in  Figure  12.  The  final  absolute  and  weighted  RMS  residuals  for  the  CC  times  are 
6  ms  and  4  ms.  The  Vp  and  Vs  models  are  stored  on  a  grid  that  has  a  horizontal  grid  spacing  of 
0.2°  in  latitude  and  longitude  and  a  vertical  grid  spacing  of  5  km  in  depth. 
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Figure  13.  Inversion  grid  nodes  for  P  wave  at  depths  7.5  km  to  12.5  km  for  the  adaptive-grid 
tomography. 

Earthquake  locations 

The  seismotectonic  structures  and  seismicity  in  the  Yunnan  region  are  controlled  by  the  collision 
between  the  Indian  and  Eurasia  plates.  As  shown  in  Figure  12,  there  are  a  number  of  active  faults 
in  the  region.  The  major  active  faults  include  the  Red  River  Fault  (F2),  the  Xiaojiang  Faut  (FI), 
and  the  Qiaohou-Weishan  Fault  (F7).  The  slip  rate  on  the  Red  River  fault  is  as  large  as  about  7  - 
9  mm/year  on  the  southern  segment  of  the  fault  (Su  and  Qin,  2001).  Several  large  earthquakes 
(e.g.  the  1925  M  7  Dali  earthquake)  have  occurred  along  the  Red  River  Fault.  Along  the 
Xiaojiang  Faut,  the  slip  rate  is  about  6.4  -  8.8  mm/year  (Fi,  1993;  Su  and  Qin,  2001).  The  strong 
earthquakes  associated  with  the  fault  include  the  1833  M8  Songming  earthquake  and  the  1966 
M6.5  Dongchuan  earthquake.  These  two  faults  along  with  the  Xianshuihe  Fault  and  Jinshajiang 
Fault  (in  the  north  of  the  region)  form  a  seismically-active  zone,  which  surrounds  a  rhomboid¬ 
like  block,  known  as  the  Sichuan- Yunnan  block.  In  addition  to  this  seismically-active  zone,  there 
are  several  other  active  areas  in  the  region.  Two  large  earthquakes  with  Ms  7.3  and  Ms  7.4 
occurred  in  the  Tengchong-Fongling  seismic  belt  (Rl)  on  May  29,  1976  only  about  100  minutes 
apart.  Similarly,  another  pair  of  earthquakes  of  Ms  7.4  and  Ms  7.2  took  place  within  about  12 
minutes  along  the  Fancang-Gengma  seismic  belt  (R2)  on  November  6,  1988.  An  M  7.0 
earthquake  ruptured  the  Fijiang  (R3)  seismic  belt  on  February  3,  1996.  Parallel  to  R3,  the 
Yongshen-Ninglang-Muli-Juilong  seismic  belt  (R4)  is  related  to  a  historical  event  of  M  7.8  in 
1515.  There  have  been  a  number  of  large  earthquakes  (e.g.  the  1970  Ms  7.8  Tonghai  earthquake) 
that  have  occurred  along  the  Tonghai- Shiping  seismic  belt  (R5). 
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Figure  12  shows  the  epicenter  distribution  for  all  the  5 144  events  relocated  in  this  study  for  the 
period  from  1999  to  2004.  For  comparison,  we  zoomed  in  on  a  subregion  to  show  catalog 
locations  and  relocated  earthquake  locations  (Figure  14).  In  general,  the  seismicity  pattern 
illustrated  by  the  original  catalog  locations  is  diffuse  and  is  very  likely  biased.  The  original 
catalog  depths  for  many  events  are  set  to  zero  due  to  the  lack  of  constraints  in  the  depth 
estimation.  A  much  sharper  image  of  the  seismicity  is  obtained  after  relocation  as  shown  in 
Figure  14.  From  Figure  12  we  can  see  that  the  majority  of  events  are  distributed  on  the  five 
seismic  belts  discussed  above.  We  can  identify  many  clear  linear  features  and  clusters,  which  are 
hard  to  identify  in  the  original  locations.  Figure  15  shows  the  focal  depth  distribution  for  all  the 
relocated  events.  From  the  figure,  we  can  see  that  the  number  of  events  increases  with  depth 
from  -2  to  12  km,  peaks  around  12  km  and  then  gradually  decreases.  Some  of  the  events 
shallower  than  2-3  km  below  surface  may  well  be  quarry  blasts.  The  majority  of  events  are 
located  within  the  depth  range  from  5  to  20  km.  This  range  likely  represents  the  seismogenic 
zone  in  the  region.  Yang  et  al.  (2005)  also  reported  that  most  of  the  events  in  the  region  are 
shallower  than  20  km.  However,  there  are  many  more  events  above  5  km  in  their  results.  The 
average  crustal  thickness  in  the  region  is  about  45-50  km.  It  is  evident  that  almost  all  the  events 
in  the  region  are  crustal  events  and  most  of  them  occur  within  the  upper  crust.  In  the  following 
we  will  discuss  the  hypocentral  distribution  for  three  seismic  belts  and  an  aftershock  cluster  in 
the  region. 
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Figure  14.  Comparison  of  catalog  and  relocated  event  locations. 
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Figure  15.  Focal  depth  distribution  for  all  the  relocated  events  in  the  Yunnan  region. 

Lancang-Gengma  seismic  belt  (R2) 

Figure  16  shows  the  relocated  seismicity  along  the  Lancang-Gengma  seismic  belt.  Most  of  the 
events  are  located  on  the  southwest  side  of  the  surface  fault  traces.  The  focal  depth  distribution 
of  events  along  the  cross  section  (A-B  in  Figure  16a)  displays  a  slight  SW  dip.  It  is  consistent 
with  the  SW-trending  fault  dips  of  ~77°  -79°  revealed  by  focal  mechanisms  of  two  large 
earthquakes  (Yang  et  al  2005,  and  references  therein).  The  focal  depths  are  centered  in  the  10-15 
km  depth  range.  There  are  only  a  few  events  deeper  than  25  km  in  the  area. 

Lijiang  seismic  belt  (R3)  and  Yongshen-Ninglang-Muli-Juilong  seismic  belt  (R4) 

The  Lijiang  and  Yongshen-Ninglang-Muli-Juilong  seismic  belts  are  two  parallel  active  seismic 
belts  within  the  Sichuan  -Yunnan  rhomboid-shaped  tectonic  block.  The  Lijiang  seismic  belt  was 
previously  described  to  be  N-S  trending  and  bounded  by  the  Lijiang  basin  (Yang  et  al.,  2005). 
More  events  along  the  northeast  part  of  the  fault  zones  are  relocated  in  this  study,  as  shown  in 
Figure  17a.  This  belt  could  be  extended  further  towards  the  northeast  along  the  strikes  of  the 
fault  traces.  Parallel  to  R3,  the  Yongshen-Ninglang-Muli-Juilong  seismic  belt  runs  northeast 
from  26.2°  to  27.2°N  and  beyond.  The  northern  segment  of  the  belt  is  out  of  the  study  region. 
These  two  NE  running  belts  are  consistent  with  the  strikes  of  major  faults  in  the  area.  There  are 
no  surface  fault  traces  around  the  segment  of  the  belt  above  26.8°N.  The  seismicity  in  this  part  of 
the  belt  likely  illustrates  a  blind  active  fault.  We  can  find  a  clear  seismicity  gap  between  these 
two  belts.  From  Figure  17b,  the  gap  can  also  be  found  in  the  focal  depth  distribution  along  the 
cross  section  (A-B).  The  majority  of  the  events  in  these  two  belts  are  within  the  focal  depth  range 
of  5-20  km. 
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Figure  16.  Hypocentral  distribution  of  relocated  seismicity  (1999-2004)  along  the  Lancang- 
Gengma  seismic  belt,  (a)  Epicentral  distribution  of  the  relocated  events,  (b)  Focal  depth 
distribution  along  the  vertical  cross  section  indicated  by  A-B  in  (a). 


20 


Figure  17.  Hypocentral  distribution  of  relocated  seismicity  along  two  parallel  belts,  Lijiang  seismic 
belt  (R3)  and  Yongshen-Ninglang-Muli-Juilong  seismic  belt,  (a)  Epicentral  distribution  of  relocated 
events,  (b)  Focal  depth  distribution  along  the  vertical  cross  section  indicated  in  (a)  by  A-B. 
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An  aftershock  cluster  adjacent  to  R3 

We  can  identify  many  earthquake  clusters  from  the  relocated  seismicity.  One  of  the  large  clusters 
is  located  just  south  of  the  Lijiang  seismic  belt.  There  are  more  than  160  events  distributed  in 
about  a  5km  x  10km  area  and  extending  from  8km  to  24  km  in  depth  (Figure  18).  Since  these 
events  occurred  within  several  days,  they  are  likely  the  aftershocks  of  an  M  6  earthquake  that 
took  place  on  November  22,  2003.  The  narrowly  distributed  event  locations  for  the  cluster 
indicate  that  the  accuracy  for  relative  locations  in  this  study  could  be  much  less  than  a  few 
kilometers. 
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Figure  18.  Relocated  hypocentral  distribution  of  an  earthquake  cluster,  (a)  Epicentral 
distribution  of  relocated  events  within  the  cluster  (b)  Focal  depth  distribution  along  the  vertical 
cross  section  indicated  in  (a)  by  A-B. 
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Vp  and  Vs  models 


Figure  19  shows  horizontal  slices  of  the  3D  P-wave  velocity  model  at  depths  of  0,  5,  10,  15,  20, 
25,  30,  40  and  45  km.  These  horizontal  slices  show  strong  lateral  velocity  heterogeneities, 
consistent  with  the  nature  of  active  fault  zones  in  this  region.  Strong  velocity  contrasts  are 
evident  across  major  faults,  such  as  the  Xiaojiang  Fault  and  Red  River  Fault,  as  shown  in  the  W- 
E  cross-sections  of  the  Vp  model  (Figure  20).  Similar  features  can  also  be  found  in  the  horizontal 
slices  and  cross-sections  of  the  3D  S-wave  model  shown  in  Figures  21  and  22. 

At  shallow  depths  (0-5km),  the  distribution  of  low  and  high  velocity  zones  is  closely  related  to 
variations  in  thickness  of  the  sedimentary  layer  and  the  regional  geology.  There  is  a  low  velocity 
anomaly  zone  to  the  west  of  longitude  99°,  along  the  Zhongdian-Lijiang-Yunlong-Wanding  zone 
(Figures  19  and  21).  Deep  Seismic  Sounding  (DSS)  studies  showed  about  3  km  of  sediments 
along  this  zone  (Hu  et  ah,  1986;  Lin  et  ah,  1993).  Between  the  Tonghai-Chuxiong  Fault  (F8)  and 
the  Nujiang  Fault  (F4),  there  is  a  large  low  velocity  anomaly  where  sediments  are  1-3  km  thick 
(He  et  ah,  2005).  This  region  is  located  around  the  junction  of  the  Yangtze  platform,  the  South 
China  fold  system  and  the  Sanjiang  fold  system  and  large  tectonic  movement  is  expected.  To  the 
east  of  the  Lancangjiang  Fault  (F4)  and  between  latitude  25°  and  28°N,  there  is  a  large  high 
velocity  anomaly  zone  where  the  sedimentary  layer  is  very  thin  (He  et  ah,  2005).  At  the  depth  of 
10  km,  the  high  and  low  velocity  anomaly  zones  are  similar  to  those  shown  at  shallow  depths. 
The  shape  of  the  high  velocity  zones  changes  very  little  from  the  surface  to  the  depth  of  10  km. 

Starting  from  the  depth  of  15  km,  we  see  low  velocity  anomalies  associated  with  the  major  fault 
zones,  such  as  the  Xiaojiang  Fault  (FI),  the  Tonghai-Chuxiong  Fault  (F8),  and  the  Nandinghe 
Fault  (F5).  At  shallower  depths,  the  velocities  around  these  zones  are  relatively  higher.  On  the 
other  hand,  the  velocity  is  higher  along  the  Zhongdian-Lijiang-Yunlong-Wanding  zone, 
indicating  the  basement  layer  velocity  in  this  zone  is  higher  than  surrounding  areas,  as  found 
from  the  phase  velocity  map  (He  et  ah,  2005).  The  low-velocity  zone  between  the  Wuliangshan 
Fault  (F3)  and  the  Red  River  Fault  (F2)  extends  all  the  way  down  to  a  depth  of  around  50  km. 

In  the  cross-sections  shown  in  Figures  20  and  22,  we  can  see  that  most  of  the  faults  are 
associated  with  a  velocity  contrast.  In  the  cross-section  along  latitude  27°N,  to  the  west  of  the 
Nujiang  Fault  (F6),  the  velocities  are  lower.  Between  the  Nujiang  Fault  (F6)  and  the  Zemuhe 
Fault,  there  is  a  broad  zone  of  low  velocity  that  is  the  central  part  of  the  Sichuan- Yunnan 
rhomboid  block.  Most  of  the  earthquakes  are  concentrated  around  the  edge  of  the  velocity 
contrast.  To  the  east  of  the  Anninghe  Fault  (F9),  there  is  a  high  velocity  anomaly  around  latitude 
27°N  starting  from  a  depth  of  25  km.  This  high  velocity  body  likely  depicts  the  Kangdian  uplift, 
related  to  an  intrusion  of  deep  material.  In  the  cross-section  along  latitude  27°N,  there  is  a  low 
velocity  depression  between  the  Lancangjiang  Fault  (F4)  and  the  Red  River  Fault  (F2).  There  is  a 
cluster  of  earthquakes  around  longitude  101°E  and  there  is  no  fault  trace  on  the  surface  around 
this  cluster.  However,  we  see  these  earthquakes  are  located  along  a  velocity  contrast  that  is 
higher  to  the  east.  It  is  likely  this  cluster  of  earthquakes  is  associated  with  a  blind  fault.  Along  the 
cross-section  at  latitude  25°N,  from  west  to  east,  it  crosses  the  Sanjiang  fold  zone,  the  Dianzhong 
depression  and  the  Kangdian  uplift.  Clear  velocity  contrasts  are  associated  with  the 
Lancangjiang  Fault,  the  Red  River  Fault  and  the  Xiaojiang  Fault.  Tengchong  Volcano  is  located 
around  longitude  98.3°  and  to  its  west  down  to  the  depth  of  50  km,  there  is  a  low  velocity 
anomaly. 
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Figure  19.  Horizontal  slices  of  the  three-dimensional  Vp  model  at  depths  of  0,  5,  10,  15,  20,  25, 
30,  35,  40,  45  and  50  km  (continued  on  the  next  page). 
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Figure  19  (continued). 
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Figure  20.  W-E  cross-sections  of  the  Vp  model  at  latitudes  of  27°,  26°  and  25°N. 
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Figure  21.  Horizontal  slices  of  the  three-dimensional  Vs  model  at  depths  of  0,  5, 10,  15,  20,  25, 
30,  35,  40,  45  and  50  km  (continued  on  the  next  page). 
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Figure  21  (continued). 
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Figure  22.  W-E  cross-sections  of  Vs  model  at  latitudes  of  27°,  26°  and  25°N. 
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7.  Earthquake  locations  and  3D  velocity  models  in  the  Sichuan  region 

For  the  Sichuan  region,  we  collected  catalog  picks  for  -5000  earthquakes  observed  on 
-50  stations.  Figure  23  shows  the  event  and  station  distribution  in  this  study.  Note  that  the 
earthquake  locations  are  the  final  product  of  our  seismic  tomography  study,  not  the  original 
catalog  locations.  There  are  -39600  P  picks  and  -37100  S  picks  in  our  data  set.  Figure  24  shows 
the  ray  path  coverage  for  P  and  S  waves.  For  the  selected  events,  there  are  at  least  8  observations 
including  P  and  S  times.  From  these  absolute  times,  we  constructed  280,000  P  and  270,000  S 
differential  times.  The  number  of  average  links  per  event  pair  is  1 1.  The  average  hypocentral 
separation  for  the  linked  event  pairs  is  -8.4  km.  Though  we  also  calculated  cross-correlation 
times  for  -300  events,  they  have  minimal  effect  on  final  event  locations  and  velocity  models 
because  of  the  limited  number  of  observations. 

We  first  solved  for  a  ID  model  at  depths  from  0  to  80  km  with  an  interval  of  5  km  using 
the  absolute  P  and  S  times.  Then  we  used  this  ID  model  as  our  initial  model  to  simultaneously 
determine  event  locations  and  velocity  models  using  both  absolute  and  differential  P  and  S  times. 
The  horizontal  grid  spacing  is  0.5°  in  latitude  and  longitude.  The  absolute  and  weighted  arrival 
time  RMS  residuals  at  the  beginning  of  the  inversion  are  1.877  s  and  4.786  s  and  they  decrease  to 
325  ms  and  53  ms  at  the  end  of  the  inversion,  respectively.  Very  large  smoothing  weight  of  600 
is  applied  because  the  data  are  noisy  and  the  accuracy  of  arrival  times  is  not  good  (-0. 1  s). 


Figure  23.  Event  and  station  distribution  for  the  Sichuan  region.  For  fault  names,  refer  to  Figure  1. 
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Figure  24.  Ray  path  coverage  for  (left)  P  and  (right)  S  waves.  Earthquakes  are  indicated  by  red 
dots  and  stations  are  indicated  by  green  circles. 

Earthquake  locations 

(1)  Xianshuihe  seismic  belt  (FI 2) 

The  Xianshuihe  Fault  is  the  northeast  boundary  of  the  Sichuan- Yunnan  rhomboid  tectonic  block. 
The  slip  rate  along  the  fault  is  about  8-15  mm/yr  (Su  and  Qin,  2001).  Although  a  series  of  large 
earthquakes  (e.g.  the  1973  M  7.6  earthquake  and  the  1981  Ms  7.0  earthquake)  occurred  along  the 
fault,  the  recent  seismicity  in  the  belt  is  relative  low  as  shown  in  Figures  23  and  25.  The 
epicentral  distribution  of  events  in  this  belt  follows  the  NW  strike  of  the  fault  trace  well.  As 
shown  in  Figure  F25b,  the  events  in  the  belt  are  usually  shallower  than  25  km. 

(2)  Longmenshan  seismic  belt  (FI  1) 

The  Longmenshan  seismic  belt  is  the  most  active  seismic  area  in  the  Sichuan  region.  As  shown 
in  Figures  23  and  26,  there  are  many  small  to  moderate  earthquakes  occurring  in  this  ~100  km 
wide  belt.  The  band  clearly  coincides  with  the  NE  strike  of  the  Longmenshan  Fault.  No 
earthquakes  with  M  >  7  have  been  documented  in  the  belt.  However,  several  M  >  6  earthquakes 
(e.g.  the  1657  M  6.5  Wenchuan  Earthquake  and  the  1970  M  6.2  Dayi  Earthquake)  occurred 
along  the  fault  historically.  The  focal  depth  distribution  along  the  cross  section  A-B  in  Figure 
26a  is  shown  in  Figure  26b.  The  event  depths  are  centered  on  ~  20  km  and  extend  to  ~40  km 
deep.  Apparently,  there  are  more  deep  events  along  the  belt  than  the  seismic  belts  in  the  Yunnan 
region. 
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Figure  25.  Relocated  seismicity  along  the  Xianshuihe  seismic  fault  (F3).  (a)  Epicentral 
distribution  of  relocated  events,  (b)  The  focal  depth  distribution  of  relocated  events  along  the 
cross  section  A-B. 
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Figure  26.  Relocated  seismicity  along  the  Longmenshan  seismic  belt  (F5).  (a)  Epicentral 
distribution  of  relocated  events,  (b)  The  focal  depth  distribution  of  relocated  events  along  the 
cross  section  A-B. 
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Vp  and  Vs  models 

Figures  27  and  28  show  horizontal  slices  of  the  3D  Vp  and  Vs  models  at  different  depths, 
respectively.  The  regions  that  are  not  well  resolved  are  shown  as  blank.  Strong  velocity 
variations  can  be  seen  in  both  the  Vp  and  Vs  models  indicating  significant  structure 
heterogeneities  in  this  region.  At  shallow  depths  of  0  and  5  km,  we  see  clear  low  velocity 
anomalies  associated  with  the  Sichuan  Basin  to  the  east,  the  northern  end  of  Xianshuihe  Fault  to 
the  west,  and  the  southern  end  of  the  Xianshuihe  Fault  and  the  northern  end  of  the  Anninghe 
Fault.  There  are  few  earthquakes  occurring  in  these  low  velocity  anomaly  zones.  The 
Longmenshan  Fault  separates  a  positive  anomaly  to  the  west  and  a  negative  anomaly  to  the  east. 
A  high  velocity  anomaly  corresponds  to  Minshan-Longmenshan  upwarping  and  Maerkang 
upwarping  (Xu  et  al.,  1985).  The  Longmenshan  Fault  is  considered  as  a  tectonic  boundary 
between  the  Songpan-Garze  Fold  system  to  the  west  and  Yangtze  Craton  to  the  east.  The 
Songpan-Garze  fold  system  is  composed  of  slightly  metamorphosed  rocks  thrust  over 
sedimentary  formation  in  the  Yangtze  Platform  (Ren  et  al.,  1980;  Wang  et  al.,  2003).  There  is 
also  an  evident  high  velocity  anomaly  located  around  latitude  30°  and  longitude  102°  from 
Panzhihua  to  Xiangcheng.  This  high  velocity  anomaly  likely  corresponds  to  a  passively 
reactivated  ancient  rift  zone  (Teng,  1987;  Huang  et  al.,  2002),  which  contains  mainly 
metamorphic  complex  and  basic  and  ultra-basic  rocks  that  are  considered  to  consist  of  high 
density  upper  mantle  materials  (Xu  et  al.,  1985). 

At  depths  of  10  to  20  km,  strong  low  velocity  anomalies  are  still  associated  with  the  Xianshuihe 
Fault  and  the  northern  segment  of  the  Anninghe  Fault.  However,  the  Longmenshan  Fault  zone  is 
associated  with  high  velocity  anomalies.  Most  of  earthquakes  are  located  on  the  boundary 
between  high-V  and  low-V  zones. 

At  deeper  depths  of  25  to  35  km,  strong  low  velocity  anomalies  are  located  in  the  southern  and 
northern  parts  of  the  study  areas  corresponding  to  the  Xianshuihe  Fault  and  the  Anninghe  Fault 
zones,  but  are  separated  by  the  Longmenshan  Fault.  This  low  velocity  anomaly  is  likely  caused 
by  high  temperature  anomalies  and  fluid  reservoirs  (Huang  et  al.,  2002)  and  may  be  weakened 
by  the  uplifting  and  folding  due  to  Indo-China  movement.  On  the  other  hand,  the  velocities  are 
higher  underneath  the  Yangtze  Platform  that  is  more  stable. 

The  low  velocity  zones  associated  with  the  southern  segment  of  the  Xianshuihe  Fault  and  the 
northern  segment  of  the  Anninghe  Fault  continue  to  a  depth  of  45  km.  In  fact,  the  low  velocity 
anomaly  still  exists  at  least  down  to  the  depth  of  65  km.  This  result  is  consistent  with  the  low  Vp 
anomaly  beneath  southeast  Tibet  from  the  teleseismic  P-wave  travel  time  tomography  (Li  et  al., 
2006)  and  the  low  phase  velocity  or  S-wave  velocity  beneath  southeast  Tibet  from  ambient  noisy 
tomography  (Yao  et  al.,  2006).  This  low  velocity  anomaly  may  be  related  to  lower  or  middle 
crustal  flow  in  the  geodynamical  models  by  Royden  et  al.  (1997)  and  Beaumont  et  al.  (2004). 
These  low  velocity  zones  may  cause  the  weakening  of  the  seismogeneic  layer  in  the  upper  crust 
and  thus  the  large  earthquakes  are  more  likely  to  happen  due  to  the  effect  of  the  tectonic  stress 
(Zhao  et  al.,  2000). 


32 


Latitude  Latitude 


depth=5  km 


depth=0  km 


1(  


101  102  103  104  105 
Longitude 


32 


31 

<L> 

=>  30 


t  .  V 


• 


\ 


■3 


*V, 

*  ‘  J 

"  n 


(O 


Sr’S. 

p$*'  ■! 


4.6  4  8  5  5.2  5.4 


29 

28 


*y 


27 


^  It;  ti  ' 

■  \hr  ' 

A 


:v 


0 

31  "21L 


101  102  103  104  105 
Longitude 


4.8  5  5.2  5.4  5.6 


32  v 
31 
30 


depth=10  km 

r3~TF2 


VJ  . 


•  •  v-.v*  ^  -  /•<** 

\  '--'M 

if  /  ?. 


r  . 


24*'^7»'v 


'101  102  103  104  105 
Longitude 


5  5.2  5.4  5.6  5.8 


depth=15  km 


1 — •{ 


§  30  ’ 

^  29  r 

'4-  >L i  i  (Jw  / 

28  ^  ft  ’#  •  • 


28  kite? 

r-  " 1  £?■  >: 

'101  102  103  104  105 
Longitude 


5.4  5.6  5.8 


Figure  27.  Horizontal  slices  of  the  3D  Vp  model  at  depths  from  0  to  75  km  at  an  interval  of  5  km 
(continued). 
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Figure  27  (continued). 
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Figure  27  (continued). 
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Figure  27  (continued). 
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Figure  28.  Horizontal  slices  of  the  3D  Vs  model  at  depths  from  0  to  45  km  at  an  interval  of  5  km 
(continued). 
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Figure  28  (continued). 
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Figure  28  (continued). 


8.  Seismic  attenuation  models  for  the  Yunnan  region 

Generally,  the  attenuation  tomography  system  is  solved  using  an  iterative  method  in  which  the  Q 
model  is  updated  according  to  t*  residuals  from  an  initial  Q  model.  We  found  out  the  Q  values 
could  be  negative  and  they  also  depend  on  the  starting  Q  values  to  some  extent.  Once  the  event 
locations  and  the  velocity  model  are  known,  the  relationship  between  t*  values  and  the  reciprocal 
Q  structure  (1/Q)  is  linear  as  shown  in  Equation  (2).  For  this  reason,  we  modified  the  program  by 
using  a  nonnegative  least  squares  (NNLS)  algorithm  (Lawson  and  Hanson,  1974)  to  directly 
solve  the  system  for  1/Q  without  iterations.  This  method  guarantees  that  the  Q  values  are 
nonnegative  and  not  dependent  on  the  initial  values.  As  noticed  by  Roth  and  Wiens  (1999),  the 
NNLS  algorithm  may  yield  some  very  high  Q  values  but  the  nonnegative  constraint  removes  the 
negative  models  that  are  not  realistic. 

We  first  use  the  -10200  P-wave  t*  values  estimated  from  an  Q-square  source  model  from  -3000 
events  on  26  stations.  The  horizontal  and  vertical  grid  intervals  are  0.5°  and  10  km.  The  velocity 
model  and  event  locations  are  taken  from  the  previous  velocity  tomography  results.  Similar  to 
the  velocity  tomography,  smoothing  and  damping  are  applied  to  regularize  the  inversion  system 
but  the  smoothing  is  applied  to  neighboring  1/Q  values,  not  perturbations.  Smoothing  weight  and 
damping  values  are  selected  based  on  trade-off  curves  between  data  variance  and  model  variance 
The  RMS  t*  value  is  0.0294  and  after  the  inversion  the  RMS  t*  residual  is  0.0101,  indicating  the 
t*  values  are  successfully  fit  to  about  2/3  of  their  true  values  with  the  resulting  Q  models. 

Figure  29  shows  horizontal  slices  of  the  3D  Qp  model  at  depths  of  0,  10,  20,  and  30  km.  At  the 
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depth  of  0  km,  the  Qp  values  are  generally  low  indicating  relatively  high  attenuation  because  of 
sediments  near  the  surface.  At  the  depth  of  10  km,  we  can  see  high  Qp  values  correlate  relatively 
well  with  high  velocity  anomalies  where  the  sedimentary  layer  is  thin,  as  shown  in  Figure  19. 
There  is  a  low-Qp  zone  between  the  Red  River  fault  (F2)  and  the  Wuliangshan  fault  (F3), 
extending  to  a  depth  of  at  least  20  km,  consistent  with  thick  sedimentary  layers  in  the  region.  The 
Kangdian  uplift  to  the  east  of  the  Anninghe  Fault  (F9)  also  shows  high  Qp  anomalies  related  to 
an  intrusion  of  deep  material. 

We  then  updated  this  Qp  model  using  differential  t*  values  measured  from  a  least  squares  line 
fitting  to  spectral  ratios  for  two  events  and  two  stations.  There  are  -127000  differential  t*  values 
corresponding  to  -1550  events  and  23  stations.  For  the  triple  difference  attenuation  tomography 
method,  only  regions  near  the  sources  and  stations  are  resolved  because  rays  for  two  events  and 
two  stations  overlap  in  the  regions  away  from  events  and  stations.  For  this  reason,  we  used  the 
Qp  model  from  conventional  attenuation  tomography  as  an  a  priori  model  in  the  triple  difference 
attenuation  tomography.  Figure  30  shows  horizontal  slice  of  the  new  Qp  model  at  depths  of  0,  10, 
20  and  30  km.  Compared  to  the  Qp  model  shown  in  Figure  29,  there  are  some  obvious  changes 
in  some  regions  where  earthquakes  are  concentrated.  For  example,  at  the  depth  of  10  km,  the 
event  cluster  around  latitude  26°N  and  longitude  100°E  was  originally  located  in  a  high  Qp 
region  (Figure  29),  but  now  it  is  more  closely  related  to  a  moderate  and  low  Qp  region  (Figure 
30).  For  comparison,  we  also  calculated  the  Qp  model  using  the  differential  t*  values  obtained 
directly  from  absolute  t*  values.  As  we  noticed  previously  that  the  two  sets  of  differential  t* 
values  correlate  very  well  (Figure  1 1),  the  Qp  models  are  very  similar  too  (Figures  30  and  3 1). 

For  the  Qs  attenuation  tomography,  we  measured  -8860  absolute  S-wave  t*  values  for  -2300 
event  on  25  stations.  The  grid  setting  is  the  same  as  for  the  Qp  attenuation  tomography.  The 
RMS  t*  value  is  0.0286  and  the  RMS  residual  after  the  inversion  is  0.0096,  indicating  the  t* 
values  are  successfully  fit  to  about  2/3  of  their  true  values  with  the  resulting  Qs  model.  Figure  32 
shows  horizontal  slices  of  the  Qs  model  at  depths  of  0,  10,  20  and  30  km.  The  high  and  low  Q 
anomalies  are  very  similar  in  Qp  and  Qs  models.  For  example,  the  high  Qs  region  around  latitude 
26°N  and  longitude  100°E  also  corresponds  to  a  high  Qp  region.  Between  the  Wuliangshan  Fault 
(F3)  and  the  Red  River  Fault  (F2),  there  also  exists  a  low  Qs  anomaly. 

We  also  updated  the  Qs  model  using  the  differential  S-wave  t*  values  measured  from  fitting 
spectral  ratios  for  two  events  and  two  stations  (Figure  33)  and  calculated  directly  from  the 
absolute  S-wave  t*  values  (Figure  34).  There  are  -122500  differential  t*  values  corresponding  to 
-1450  events  and  24  stations.  There  are  some  clear  changes  in  the  Qs  model  using  the 
differential  t*  values.  For  example,  at  the  depth  of  20  km,  in  the  western  side  of  the  study  region, 
the  Qs  values  are  uniformly  high  from  the  inversion  of  absolute  t*  values  (Figure  32).  However, 
locally  low  Qs  anomaly  zones  are  imaged  around  latitude  25°N  and  longitude  100°E  where 
earthquakes  are  distributed  (Figures  33  and  34).  Similar  features  can  be  observed  in  the  two  Qs 
models  updated  using  different  sets  of  differential  S-wave  t*  values.  Because  of  strong 
correlations  in  differential  t*  values  measured  using  the  two  different  methods  and  the 
similarities  in  the  Q  models  for  both  P  and  S  waves,  it  is  therefore  reasonable  to  use  the 
differential  t*  values  calculated  directly  from  the  absolute  t*  values. 


40 


Latitude  Latitude 


depth=0  km 


96  98  100  102  104 

depth=20  km 


96  98  100  102  104 

Longitude 


depth=10  km 


96  98  100  102  104 

depth=30  km 


_ _ _ LAu _ _ 

96  98  100  102  104 

Longitude 


Figure  29.  Horizontal  slices  of  the  Qp  model  at  depths  of  0,  10,  20,  and  30  km  calculated  using 
the  t*  values  estimated  from  the  omega-square  source  model. 
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Figure  30.  Horizontal  slices  of  the  Qp  model  at  depths  of  0,  10,  20,  and  30  km  calculated  using 
the  differential  t*  values  estimated  from  the  waveform  fitting. 
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Figure  31.  Horizontal  slices  of  the  Qp  model  at  depths  of  0,  10,  20,  and  30  km  calculated  using 
the  differential  t*  values  estimated  directly  from  absolute  t*  values. 
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Figure  32.  Horizontal  slices  of  the  Qs  model  at  depths  of  0,  10,  20,  and  30  km  calculated  using 
the  t*  values  estimated  from  the  omega-square  source  model. 
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Figure  33.  Horizontal  slices  of  the  Qs  model  at  depths  of  0,  10,  20,  and  30  km  calculated  using 
the  differential  t*  values  estimated  from  the  waveform  fitting. 
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Figure  34.  Horizontal  slices  of  the  Qs  model  at  depths  of  0,  10,  20,  and  30  km  calculated  using 
the  differential  t*  values  estimated  directly  from  the  absolute  t*  values. 
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9.  Assembling  ground  truth  events  in  the  Yunnan  region 

Determining  accurate  seismic  locations  and  estimating  their  uncertainties  are  of  fundamental 
importance  to  ground-truth  based  nuclear  explosion  monitoring.  The  traditional  ways  to 
determine  location  uncertainties  include  the  methods  based  on  a  posteriori  residual  distribution 
(Flinn,  1965)  or  a  priori  uncertainty  for  phase  picking  and  travel  time  prediction  (Evemden, 
1969).  Both  methods  assume  that  the  error  processes  are  Gaussian,  zero  mean  and  uncorrelated. 
In  most  seismic  location  practices,  these  basic  statistical  assumptions  are  not  satisfied.  Thus,  the 
uncertainty  in  these  formal  calculations  could  be  underestimated  (Myers  and  Schultz,  2001).  An 
alternative  way  to  assess  catalogue  location  accuracy  is  based  on  the  station  geometry.  Using 
Monte  Carlo  simulations  of  network  geometry,  Bondar  et  al.  (2004)  found  the  local  network 
locations  (0°-  2.5°  epicentral  distance)  are  accurate  to  within  5  km  with  a  95  percent  confidence 
level  (GT5  95  percent)  when  the  network  meets  the  following  criteria:  (1)  there  are  at  least  10 
stations  within  250  km;  (2)  an  azimuth  gap  of  less  than  1 10°;  (3)  a  secondary  azimuth  gap  of  less 
than  160°;  (4)  at  least  one  station  within  30  km  from  the  event  epicenter. 

We  use  similar  criteria  as  defined  by  Bondar  et  al.  (2004)  to  classify  the  earthquakes  with 
magnitude  4  or  greater  in  the  Yunnan  region  from  1999  to  2004.  The  Yunnan  seismic  network 
consists  of  26  stations  distributed  in  a  6°x7°  area.  Some  stations  are  just  beyond  the  local 
epicentral  distance  of  280  km.  Thus,  we  have  to  tailor  the  first  two  sets  of  criteria  to  satisfy  our 
situation.  There  are  a  number  of  events  with  magnitude  4  or  above  that  meets  the  first  3  criteria 
of  GT5  95  percent  defined  by  Bondar  et  al.  (2004).  However,  only  a  few  of  them  conform  to  the  4th 
criterion.  As  indicated  in  Bondar  et  al.  (2004),  this  constraint  is  mainly  used  to  ensure  that  the 
event  is  within  the  crust.  As  discussed  in  previous  sections,  our  candidate  events  are  all  crustal 
events  with  high  confidence,  thus,  we  reasonably  relax  the  4th  constraint  to  that  there  is  at  least 
one  station  within  100  km  of  the  epicenter.  There  are  a  total  of  22  earthquakes  with  magnitude  4 
or  greater  to  match  the  4  criteria.  The  average  number  of  stations  that  record  these  events  is  19 
with  an  average  epicentral  distance  less  than  280  km.  Taking  into  account  that  the  accurate  3D 
velocity  model  used  in  relocation,  the  location  accuracy  for  these  earthquakes  should  be  well 
within  5  km.  We  ascribe  them  to  the  GT5  category.  The  catalog  information  for  these  GT5 
reference  events  is  listed  in  Table  1  and  their  locations  are  shown  in  Figure  35.  From  Figure  35, 
we  can  see  that  these  GT5  reference  events  are  located  in  the  central  area  of  the  study  region  and 
are  surrounded  by  the  network  stations. 
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Table  1.  Event  information  for  the  assembled  GT5  earthquakes 


Year 

mon 

day 

hour 

min 

sec 

lat. 

Lon. 

depth 

Mag 

1999 

7 

17 

15 

55 

47.23 

24.7545 

99.6991 

10.30 

4.2 

1999 

12 

7 

11 

22 

42.36 

24.6331 

102.8813 

12.11 

4.0 

2000 

1 

15 

0 

20 

15.15 

25.5782 

101.0424 

9.62 

4.7 

2000 

1 

15 

1 

3 

43.28 

25.5885 

101.0179 

17.36 

4.1 

2000 

1 

23 

0 

21 

19.65 

25.5785 

101.054 

14.04 

4.3 

2000 

6 

30 

0 

27 

45.79 

25.5847 

101.0459 

13.26 

4.1 

2000 

8 

21 

13 

25 

38.65 

25.7424 

102.2824 

4.84 

5.1 

2001 

3 

13 

11 

42 

35.67 

25.0082 

101.4257 

11.20 

4.0 

2001 

4 

15 

2 

11 

57.52 

24.0376 

99.4911 

29.42 

4.4 

2001 

5 

24 

2 

23 

37.39 

24.9524 

103.0299 

18.92 

4.2 

2001 

9 

4 

4 

5 

54.47 

23.6743 

100.6269 

13.22 

4.9 

2001 

9 

6 

12 

13 

32.57 

26.1814 

100.7509 

8.21 

4.2 

2002 

5 

13 

1 

52 

3.17 

25.3018 

103.1218 

8.71 

4.6 

2002 

7 

25 

0 

34 

56.31 

25.5838 

101.0258 

14.85 

4.0 

2003 

7 

21 

15 

16 

28.94 

25.9857 

101.2480 

3.22 

6.5 

2003 

8 

12 

15 

59 

39.97 

26.0222 

101.0420 

6.70 

4.7 

2003 

8 

18 

10 

0 

16.09 

25.9577 

101.2380 

12.79 

4.7 

2003 

9 

27 

1 

6 

28.55 

25.9644 

101.2311 

15.59 

4.2 

2004 

2 

9 

4 

30 

11.40 

25.7581 

99.9375 

17.77 

4.6 

2004 

5 

14 

1 

54 

51.74 

25.8603 

102.2713 

9.44 

4.0 

2004 

6 

10 

2 

14 

24.93 

25.5386 

101.8908 

18.22 

4.5 

98  99  100 


101  102 
longitude(°) 


103  104  105 


Figure  35.  Locations  of  GT5  reference  earthquakes. 
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10.  Conclusions 


We  have  developed  an  adaptive-mesh  "sphere-in-a-box"  version  of  the  DD  tomography  code  and 
applied  it  to  the  Yunnan-Sichuan  region  using  catalog  P  and  S  picks  and  waveform  CC  times  (for 
the  Yunnan  region).  Both  P-  and  S-wave  velocity  models  show  strong  velocity  heterogeneities, 
consistent  with  the  nature  of  this  active  transition  zone  between  the  Yangtze  platform  to  the  east 
and  the  Tibetan  plateau  to  the  west.  Strong  velocity  contrasts  are  evident  across  some  major  fault 
zones  and  the  velocity  models  are  consistent  with  the  regional  geology. 

The  accuracy  of  event  locations  is  greatly  improved  compared  to  catalog  locations  due  to  the  use 
of  more  accurate  3D  velocity  models  and  more  accurate  differential  times  (waveform  CC  times 
and/or  catalog  differential  times).  A  much  sharper  image  of  the  seismicity  is  obtained  after 
relocation  for  the  Yunnan  region  compared  to  the  diffuse  and  biased  original  catalog  locations. 
Many  linear  features  and  clusters  can  be  identified  in  the  newly  determined  event  locations.  The 
majority  of  the  events  in  the  Yunnan  region  are  located  within  the  depth  range  from  5  to  20  km 
and  peak  around  12  km.  The  previous  study  using  the  DD  location  method  also  found  most  of 
the  earthquakes  in  the  Yunnan  region  are  shallower  than  20  km,  but  they  reported  a  rather 
uniform  event  distribution  from  0  to  20  km  (Yang  et  ah,  2005).  There  are  more  deeper 
earthquakes  down  to  40  km  in  the  Sichuan  region.  Using  the  modified  Bondar  criteria  to  select 
GT  events,  we  selected  22  earthquakes  with  magnitude  4  or  above  to  be  likely  GT5  events. 

We  also  measured  absolute  and  differential  P-  and  S-wave  t*  values  for  the  Yunnan  region.  Qp 
and  Qs  models  are  calculated  using  the  absolute  and  differential  t*  values  using  the  NNLS 
algorithm  to  prevent  unphysical  negative  Q  values.  Overall,  the  Qp  and  Qs  models  using  the 
absolute  t*  values  show  many  similar  features  and  are  consistent  with  the  velocity  models  and 
the  local  geology.  The  Qp  and  Qs  models  updated  using  the  differential  t*  values  show  evident 
improvement  in  regions  where  earthquakes  are  distributed.  We  also  found  that  the  differential  t* 
values  measured  using  the  least  squares  line  fit  to  the  spectral  ratios  for  two  events  and  two 
stations  correlate  well  with  those  calculated  directly  from  absolute  t*  values.  Therefore,  it  is 
reasonable  to  use  the  absolute  t*  values  to  directly  construct  the  differential  t*  values  to  update 
the  Q  models. 
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